#ifndef VECTOR2_HPP_
#define VECTOR2_HPP_

#include <cmath>
#include <iostream>
#include <string>
#include <sstream>

//#define M_PI 3.141592653589793238462643  //swig ignore

class Vector2 {
//------------------------------------------------------------------------------
//  data members
//------------------------------------------------------------------------------
protected:

   double    fX;    // components of the vector
   double    fY;
//------------------------------------------------------------------------------
//  function members
//------------------------------------------------------------------------------
public:

   Vector2 ();
   Vector2 (double *s);  // swig ignore
   Vector2 (double x0, double y0);
   virtual ~Vector2();
                                        // ****** unary operators

   Vector2&       operator  = (Vector2 const & v);
   Vector2&       operator += (Vector2 const & v);
   Vector2&       operator -= (Vector2 const & v);
   double        operator *= (Vector2 const & v);
   Vector2&       operator *= (double s);
   Vector2&       operator /= (double s);

  Vector2 operator+(const Vector2& v) {
    return Vector2(fX + v.fX, fY + v.fY); 
  }

  Vector2 operator-(const Vector2& v) {
    return Vector2(fX + v.fX, fY + v.fY); 
  }

  double operator*(const Vector2& v) {
    return fX*v.fX + fY*v.fY; 
  }

  Vector2 operator*(double s) {
    return Vector2(fX*s, fY*s);
  }

  Vector2 operator/(double s) {
    return Vector2(fX/s, fY/s);
  }

  // ****** binary operators

   friend Vector2       operator + (const Vector2&, const Vector2&); 
   friend Vector2       operator + (const Vector2&, double  );
   friend Vector2       operator + (double  , const Vector2&);
   friend Vector2       operator - (const Vector2&, const Vector2&);
   friend Vector2       operator - (const Vector2&, double  );
   friend double       operator * (const Vector2&, const Vector2&);
   friend Vector2       operator * (const Vector2&, double  );
   friend Vector2       operator * (double  , const Vector2&);
   friend Vector2       operator / (const Vector2&, double  );
   friend double       operator ^ (const Vector2&, const Vector2&);

                                        // ****** setters
   void Set(const Vector2& v);
   void Set(double x0, double y0);
   void Set(float  x0, float  y0);  //swig ignore

                                        // ****** other member functions

   double Mod2() const { return fX*fX+fY*fY; };
   double Mod () const;

   double Px()   const { return fX; };
   double Py()   const { return fY; };
   double X ()   const { return fX; };
   double Y ()   const { return fY; };

                                        // phi() is defined in [0,TWOPI]

   double Phi           () const;
   double DeltaPhi(const Vector2& v) const;
   void     SetMagPhi(double mag, double phi);

                                        // unit vector in the direction of *this

   Vector2 Unit() const;
   Vector2 Ort () const;

                                        // projection of *this to the direction
                                        // of Vector2 vector `v'

   Vector2 Proj(const Vector2& v) const;

                                        // component of *this normal to `v'

   Vector2 Norm(const Vector2& v) const;

                                        // rotates 2-vector by phi radians
   Vector2 Rotate (double phi) const;

                                        // returns phi angle in the interval [0,2*PI)
   static double Phi_0_2pi(double x);                                                                               // returns phi angle in the interval
                                        // returns phi angle in the interval [-PI,PI)
   static double Phi_mpi_pi(double x);

   std::string ToString();
};


inline Vector2& Vector2::operator  = (Vector2 const& v) {fX  = v.fX; fY  = v.fY; return *this;}
inline Vector2& Vector2::operator += (Vector2 const& v) {fX += v.fX; fY += v.fY; return *this;}
inline Vector2& Vector2::operator -= (Vector2 const& v) {fX -= v.fX; fY -= v.fY; return *this;}
inline double   Vector2::operator *= (const Vector2& v) { return(fX*v.fX+fY*v.fY); }

inline Vector2& Vector2::operator *= (double s) { fX *=s; fY *=s; return *this; }
inline Vector2& Vector2::operator /= (double s) { fX /=s; fY /=s; return *this; }

inline Vector2  operator + (const Vector2& v1, const Vector2& v2) {
   return Vector2(v1.fX+v2.fX,v1.fY+v2.fY);
}

inline Vector2  operator + (const Vector2& v1, double bias) {
   return Vector2 (v1.fX+bias,v1.fY+bias);
}

inline Vector2  operator + (double bias, const Vector2& v1) {
   return Vector2 (v1.fX+bias,v1.fY+bias);
}

inline Vector2  operator - (const Vector2& v1, const Vector2& v2) {
   return Vector2(v1.fX-v2.fX,v1.fY-v2.fY);
}

inline Vector2  operator - (const Vector2& v1, double bias) {
   return Vector2 (v1.fX-bias,v1.fY-bias);
}

inline Vector2  operator * (const Vector2& v, double s) {
   return Vector2 (v.fX*s,v.fY*s);
}

inline Vector2    operator * (double s, const Vector2& v) {
   return Vector2 (v.fX*s,v.fY*s);
}

inline double operator * (const Vector2& v1, const Vector2& v2) {
   return  v1.fX*v2.fX+v1.fY*v2.fY;
}

inline Vector2     operator / (const Vector2& v, double s) {
   return Vector2 (v.fX/s,v.fY/s);
}

inline double   operator ^ (const Vector2& v1, const Vector2& v2) {
   return  v1.fX*v2.fY-v1.fY*v2.fX;
}

inline  double Vector2::DeltaPhi(const Vector2& v) const { return Phi_mpi_pi(Phi()-v.Phi()); }

inline  Vector2 Vector2::Ort () const { return Unit(); }

inline  Vector2 Vector2::Proj(const Vector2& v) const { return v*(((*this)*v)/v.Mod2()); }

inline  Vector2 Vector2::Norm(const Vector2& v) const {return *this-Proj(v); }

inline void Vector2::Set(const Vector2& v   )     { fX = v.fX; fY = v.fY; }
inline void Vector2::Set(double x0, double y0) { fX = x0  ; fY = y0 ;  }
inline void Vector2::Set(float  x0, float  y0)     { fX = x0  ; fY = y0 ;  }


double const  kPI        = M_PI;       //swig ignore
double const  kTWOPI     = 2.*kPI;     //swig ignore


Vector2::Vector2()
{
   //constructor
   fX = 0.;
   fY = 0.;
}


Vector2::Vector2(double *v)
{
   //constructor
   fX = v[0];
   fY = v[1];
}


Vector2::Vector2(double x0, double y0)
{
   //constructor
   fX = x0;
   fY = y0;
}


Vector2::~Vector2()
{
}


double Vector2::Mod() const
{
   // return modulo of this vector
   return sqrt(fX*fX+fY*fY);
}


Vector2 Vector2::Unit() const
{
   // return module normalized to 1
   return (Mod2()) ? *this/Mod() : Vector2(); 
}


double Vector2::Phi() const
{
   // return vector phi
   return M_PI+atan2(-fY,-fX); 
} 


double Vector2::Phi_0_2pi(double x) {
   // (static function) returns phi angle in the interval [0,2*PI)
  if(x != x) // portable way to check {
    {
     std::cerr << "Vector2::Phi_0_2pi: called with NaN" << std::endl;
     return x;
    }
  
  while (x >= kTWOPI) x -= kTWOPI;
  while (x <     0.)  x += kTWOPI;
  return x;
}


double Vector2::Phi_mpi_pi(double x) {
   // (static function) returns phi angle in the interval [-PI,PI)
   if(x != x){
     std::cerr << "Vector2::Phi_mpi_pi: called with NaN" << std::endl;
     return x;
   }
   while (x >= kPI) x -= kTWOPI;
   while (x < -kPI) x += kTWOPI;
   return x;
}


Vector2 Vector2::Rotate (double phi) const
{
   //rotation by phi
   return Vector2( fX*cos(phi)-fY*sin(phi), fX*sin(phi)+fY*cos(phi) );
}


void Vector2::SetMagPhi(double mag, double phi) 
{
   //set vector using mag and phi
   double amag = fabs(mag);
   fX = amag * cos(phi);
   fY = amag * sin(phi);
}

std::string Vector2::ToString() {
  std::stringstream ss;
  ss << "Vector2(" << fX << "," << fY << ")";
  return ss.str();
}




#endif
